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(54) Abstract Title 

Attenuating multiples in marine seismic data 

(57) Seismic wavefield data is recorded by seismic multi-component receivers, the data having 
measurements relating to pressure (p) and velocity (v , v , v ) components of the wavefield. The method 
involves choosing a component of the wavefield and selecting two or more components of the wavefield data 
necessary to attenuate multiples in the chosen component. One or more of the selected components is filtered 
using a spatial filter which combines data components measured at different receiver locations and operates 
indiscriminately on P- and S- waves. The selected and filtered selected components are then summed to 
generate a reduced multiple variant of the chosen component. 
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Multiple attenuation of mult i -component sea-bottom data 

The present invention relates to methods of reducing or 
5 attenuating multiples in multi -component sea-bottom data. 

BACKGROUND OF THE INVENTION 

10 Marine seismic surveys are usually conducted by towing an energy 
source and seismic detectors behind a vessel. The source imparts 
an acoustic wave to the water, creating a wavefield which 
travels coherently into the underlying earth. As the wavefield 
strikes interfaces between earth formations, or strata, it is 

15 reflected back through the earth and water to the detectors, 
where it is converted to electrical signals and recorded. 
Through analysis of these signals, it is possible to determine 
the shape, position and lithology of the sub-bottom formations. 
In other marine survey methods, the detectors and/or sources are 

20 placed at or close to the sea bottom or in wells. 

A problem encountered in marine surveying - as well as in 
inverse vertical seismic profiling or *VSP* - is that of water 
column reverberation. The problem, which arises as a result of 

25 the inherent reflectivity of the water surface and bottom, may 
be explained as follows. A seismic wave generated in (or 
reflected of) earth strata passes into the water in a generally 
upward direction. This wave, termed the "primary*, travels 
through the water and past the seismic detector which records 

30 its presence. The wavefield continues upward to the water's 

surface, where it is reflected back downwards. This reflected, 
or "ghost*, wavefield also travels through the water and past 
the detector(s), where it is again recorded. Depending upon the 

nature of the earth material at the water's bottom, the ghost 
35 wavefield may itself be reflected upwards through the water, 
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giving rise to a series of one or more subsequent ghost reflec- 
tions or multiples. 

Free-surface multiple reflections can be classified according to 
5 their order, which is equal to the number of reflections from 
the free-surface. Thus, first order free-surface reflections 
comprise energy initially travelling downwardly from the sources 
(as opposed to "ghosting" where energy travels upwardly and is 
reflected from the free surface) , is reflected upwardly from the 
10 sea bed or a boundary below the sea bed, and is then reflected 
downwardly from the free-surface to the hydrophones. Second 
order free-surface multiple reflections undergo two downward 
reflections from the sea-surface before being detected by the 
hydrophones, and so on. 

15 

This reverberation of the seismic wavefield in the water 
obscures seismic data, amplifying certain frequencies and 
attenuating others, thereby making it difficult to analyse the 
underlying earth formations. 

20 

In instances where the earth material at the water bottom is 
particularly hard, excess acoustic energy or noise generated by 
the seismic source can also become trapped in the water column, 
reverberating in the same manner as the reflected seismic waves 
25 themselves. This noise is often high in amplitude and, as a 
result, rends to cover the weaker seismic reflection signals 
sought for study. 

In the art, Ruehle, in U.S. Patent No. 4,486,865, discloses a 
3 0 technique for reducing ghosting wherein a pressure detector and 
a particle velocity detector positioned in close proximity to 
one another in the water. The output of at least one of the 
detectors is gain-adjusted and filtered, using a deconvolution 
operation having a predetermined amount of white noise to the 
35 zero lag of the autocorrelation function. The patent suggests 
that, by adding this deconvolved/gain-adjusted signal to the 
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output of the other detector, ghost reflections may be 
cancelled. 

Dragoset, in U.S. Patent No. 5,365,492, describes an iterative 
5 method of determining an scaling factor used to sum simultaneous 
pressure and velocity measurements. 

Haggerty, in U.S. Pat. No. 2,757,356, discloses a marine seismic 
reflection surveying system in which two seismometer spreads are 
10 disposed at two distinct depths in the water such that water 
column reverberations received by them are 180 degrees out of 
phase. By combining the output of the detectors, the patent 
suggests that the reverberations will cancel. 

15 White, in 'Seismic Wave Radiation - Transmission and 

Attenuation', McGraw-Hill, 19 65, and later Barr, in U.S. Pat. 
No. 4,979,150, both propose the use of hydrophone /geophone pair 
to separate upwardly and downwardly travelling waves. Both use a 
plane wave decomposition of the wave equation to determine a 

2 0 scaling factor which allows to add measurements from the 
different detectors types. The suggested solution assumes 
however normal incidence of the waves. 

Ikelle et al . , in the UK Patent Application GB-A-9710435 , use a 
25 Born series approach to eliminate multiples from sea-bottom 
surveys . 

Berni, in U.S. Pat. No. 4,345,473, suggests the use of a 
vertical component accelerometer in combination with 
30 hydrophone for cancelling surface-reflected noise in 
marine seismic operations. 

Gal'perin, in 'Vertical Seismic Profiling,* Special Publication 
No. 12 of the Society of Exploration Geophysicists , suggests the 
35 use of a seismic detector which combines the output of a 
pressure and velocity sensor for use in VSP surveying. 
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Moldoveanu, in U.S. Pat. No. 5,621, 700, uses a pressure and 
velocity sensitive sea-bottom cable. The reverberations are 
attenuated by an averaging process involving the steps of adding 
5 the product of the pressure data times the absolute value of the 
velocity data and the product of the velocity data times the 
absolute value of the pressure data. 

L. Amundsen and A . Reitan published in Geophysics, Vol. 60, 2, 
10 1995, 563-572 a method for decomposing multi-component sea floor 
data, employing decomposition filters determined by plane wave 
analysis . 

K.M. Schalkwijk et al . published in SEG, Expanded Abstracts, 
15 1997, p. 8-11, a method of simultaneously decomposing the 

recorded wavefield into up- and downgoing and compressional (P-) 
and shear CS-) waves in a one-step decomposition. The described 
method did not provide satisfactory results when applied to real 
data . 

20 

In view of the above-cited references, it is seen as an object 
of this invention to provide methods of eliminating ghosts or 
multiples from marine surveys, particularly sea bottom 
measurements. It is a specific object of the invention to 
25 provide such a method without making assumption about the 
direction of incidence of the acoustic waves. 



SUMMARY OF THE INVENTION 

30 

In accordance with the present invention there is provided a 
method of reducing multiples in data recorded during a marine 
seismic survey, comprising the steps of (a) providing seismic 
wavefield data recorded by seismic multi-component receivers 
35 wherein said data comprise measurements related to pressure and 
velocity components of said wavefield; (b) choosing a component 
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of the wavefield; (c) selecting two or more components of the 
provided wavefield data necessary to attenuate multiples in said 
chosen component; (d) filtering one or more of said selected 
components using a spatial filter which combines data components 
5 measured at different receiver locations and operates 

indiscriminately on P- and S-waves; and (e) suinming the selected 
and filtered selected components to generate a reduced multiple 
variant of said chosen component. 

10 In contrast to known attempts of scaling the different 

components of the wavefield before summation, such as described 
for example in the U.S. Patent Nos . 4,486,865 and 5,365,492, the 
filter of the present invention is a spatial filter using a 
convolution operation which combines several traces to generate 

15 one output trace. A trace is the measurement of one component of 
the wavefield at one location. 

In its most general embodiment, the filter is also frequency 
variant. Therefore it is possible to consider contributions of 
20 the wavefield with non-zero angle of incidence in a unified 

manner without having to separately calculate angle dependant 
scaling factors. 

Regards the method described by K.M. Schalkwijk et al . , in SEG, 
25 Expanded Abstracts, 1997, p. 8-11, the present invention 

separates the step of suppressing the downgoing multiples from 
the decomposition into compressional (P-) and shear (S-) waves. 
This approach leads to superior results when applied to real 
data . 

30 

The demultiple filter is preferably derived from a solution of 
the elastodynamic wave equation using the representation 
theorem. The demultiple filter attempts to reconstruct the 
downgoing wavefield from multi-component measurements. By 
3 5 subtracting the reconstructed downgoing wavefield from the 

measured total wavefield, the upgoing wavefield can be derived. 
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The filter is preferably a spatial filter using a convolution 
operation which combines several traces (measurements of single 
receivers) to generate one output trace. In its most general 
5 embodiment, the filter is also frequency variant. 

As the demultiple filter is derived from the representation 
theorem instead of plane wave decomposition, the sea bottom can 
be arbitrarily irregular. The filter eliminates water 
10 reverberations as the methods described in the prior art. But it 
also attenuates other multiples that are not suppressed using a 
known method. Furthermore, the new method naturally encompasses 
non-zero offset signals without having to introduce separate 
directional factor . 

15 

In practice, measurement which combine signals from two 
different receivers, such as geophones and hydrophones, require 
additional calibration to match their output and, particularly 
in case of ocean bottom measurements, to match differences in 
2 0 the coupling of receiver and sea floor. 

These and other features of the invention, preferred embodiments 
and variants thereof, and further advantages of the invention 
will become appreciated and understood by those skilled in the 
2 5 art from the detailed description and drawings below. 

BRIEF DESCRIPTION OF DRAWINGS 

30 FIG. 1 schematically illustrates marine seismic events as 

registered in a sea bottom cable; 

FIGs . 2A,B schematically illustrate up- and downgoing wavefield 

as registered in a sea bottom cable; 



35 
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FIG. 3 



graphically denotes the boundary conditions for a 
marine seismic survey using a sea bottom cable; 



5 



FIG. 4 



shows a flow diagram illustrating basic steps of a 
method in accordance with the invention. 



MODE(S) FOR CARRYING OUT THE INVENTION 



10 Referring to FIG. 1, there are shown typical sea bottom seismic 
events. A solid black circle denotes a source. A triangle 
denotes a receiver. In the examples, the sources are located 
close to the surface 10 of the sea. The receivers are located at 
the sea bottom 11. Below the sea bottom a single reflector 12 is 

15 shown at arbitrary depth. Signal paths are shown for various 
seismic events that lead to multiples in the recorded data. 
Depending on the direction of incidence at the receiver site, 
the event can be separated into upgoing (FIG.2A) and downgoing 
(FIG . 2B) events . 



As illustrated by FIG. 2A, the downgoing wavefield contains 
essentially multiples, i.e., water-layer reverberations 21 and 
receiver peg-leg multiples 22. It also includes the sea bottom 
primary 20. The upgoing wavefield, shown in FIG. 2B, contains 

2 5 all primaries 23 except the sea floor primary and some multiples 

events, in particular the source peg-leg multiples 24 which are 
not included in the downgoing wavefield. 

The new demultiple method as described below is based on the 

3 0 elastodynamic representation theorem. Let Vj be the jth component 

of the particle velocity, and p the pressure. The corresponding 
demultipled fields just below the sea floor are denoted by a - 
(tilde) superscript. It can be shown that the demultiple of Vj 
can be written as 



20 



35 
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V-j(CO, X) = Vj(CO, X) - 

[1] J dS (x* ) [ (icop (co, x' ) G 3j <x' , x, k a , k p ) - 

ViCco, x' ) n 1 (x' , x, k a , k p ) ] 

at a point x in a given domain D limited by a surface S, and 
where Gij is the Green's tensor of the elastodynamic displacement 

5 field, Dp,* is the Green's stress tensor, and ka = CD/a ) and kp = 

co/P are the P- and S-wave wavenumbers , respectively, and <B is 
the circular frequency. Using Hooke's law, a dernultipled 
pressure measurement can be written as 

p (co, x) = i(co p(x) )" 1 { (oc 2 (x) - 2 P 2 (x) ) 

10 [2] , 

[9 1 ^ 1 (co / x) + 3 2 v 2 (c0/ x) ] + or(x) d 3 v 3 (cD, x) } 

where the dernultipled velocities v are given in eq. [1] , and the 
d . denote partial derivatives, p(x) gives the density at point x. 

15 It is important to note that the equations [1] and [2] are 

general demultiple equations for ocean bottom data, and valid 
for an arbitrarily shaped sea floor and laterally varying medium 
parameters. To evaluate the equation, dynamic ray tracing can be 
used to derive a numerical expression for the Green's tensor 

20 elements. The sea floor profile S is an absorbing boundary for 
the Green's tensor. 



An analytic solution of the equations [1] and [2] can be derived 
by making simplifying assumptions. In the example below, it is 
25 for example assumed that the sea floor parameter are constant. 

For the evaluation the notation used in the equations [1] and 
[2] is modified. The points on the surface S are now denoted by 

^, and those in the volume D are denoted by x, x' . The 
30 boundaries used for evaluating the integral equation [1] are 
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shown in Figure 3, where the domain D is bounded between the sea 
bottom 31 and another reference level 32 just below the sea 
bottom. The right and left borders of domain D extend to 
infinity- By bringing the reference level very close to sea 
5 floor, the computation of the downgoing wavefield can be reduced 
to evaluating the elastodynamic representation integral along 
the sea floor. By also assuming that the medium parameters just 
below the sea floor are constant, use can be made of the free- 
space Green's tensors for Gij(x, a>, x') and Zp q)c (x, co, x' ) . With 
10 these approximations, the elastodynamic representation integral 
can be analytically evaluated. 

Let {p,v x ,v y/ v z } be the four component seismic data recorded 
directly at the sea floor; p denotes the pressure field, v x and 

15 v y denote the horizontal components of the particle 

velocity, and v 2 the vertical component of the particle velocity. 
In accordance with the standard up-down decomposition at a 
particular location, for instance the sea bottom, {p,v x ,v y ,v 2 } is 
decompose into a downgoing wavefield {p D , v x D , v y D , v 2 D } and a upgoing 

20 wavefield {p u , v x u , v y °, v 2 u } such that 

[3] {p, v x , v y , v 2 } = {p D , v x D , v y D , v 2 D }+ {p U ,v x u , v y u , v z u } 

By using the representation theorem as described earlier, it is 
25 found that the upgoing wavefield {p u , v x u , v y u , v z u } can be determined 
using 
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[4] 



(x s , co)J 



p u (x s/ co) = I [p(x s/ co) - r p;v *(x s , co) * ^ 

v x u (x s , ^ co) = I [v x (x s , co) - r v -^(x s , %. co) * v z (x s , 5, co)l 

v y u (x s , co) = I |y y (x., co) - r v - Vz (x s/ \, co) * v z (x s , \, co)] 

v z u (x s , fc, co) = I [v z (x s , \, co) - r v - v *(x s , co) * v x (x s/ \, co) - 
r Vz?Vy (x s , co) * v y (x s , ^ co) - r v ^(x s , co) * p(x s , \, co)J 



10 



where x s = (x s ,y s , z s ) is the source position and 2; = (£i, £2, £3) is the 
receiver position. The * symbol denotes spatial convolution with 
respect to The coefficients T are known analytically. For 



example : 



r p; M4, co) = -icop 



[5] 



where 



1 + 4 



r £j {di + a*) 2 g p (o, %, co) - icop 4^ J (a; + ai>| g s (o, co) 



[6] 



g P (0, t^^Sg^pfif 
g s (0, S,co) = ^|expj-i||^, 



15 and where the constants p, a and p are density, P-wave velocity 
and S-wave velocity just below the sea bottom, respectively. As 
the receivers are positioned at the sea bottom, the derivatives 
in [5] are carried out along the sea floor surface (the 

derivative variable is £) . 
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For practical applications, the filter Tmay by approximated to 
limit the filter length and hence increase the speed at which 
the method can be implemented. 

5 

Of the various methods known to reduce the length of a filter 
while preserving its efficiency/ two methods are described 
below: 

10 Using a Taylor expansion on the filter in wavenumber domain the 
filter can be written as an infinite series in the horizontal 
wavenumber. By only keeping a given number of terms, these can 
be inverse Fourier transformed to space domain giving an 
approximate spatial demultiple filter. This filter is expressed 
15 in terms of powers of second derivative operators, which again 
can be approximated to the degree of accuracy required. This 
approach is the equivalent to developing one-way wave equations 
for wavefield extrapolation/migration resulting in the so-called 
well-known 5, 15, 30, 45 etc. degree equations - handling waves 
20 propagating with up to 5, 15, 30, 45 etc. degrees from the 
vertical axis correctly. 

Using numerically optimized filters to determine by optimization 
a spatial filter of pre-defined length that has a wavenumber 
expression that fits in some sense the true wavenumber filter 
over a pre-defined wavenumber range. Such filters can be very 
short compared with the one given by eq . [4] . This approach is 
similar to the one used to optimize derivative operators for 
numerical derivation . 

Furthermore, the accuracy and the scope of the above-described 
method can be increased using various extensions: 

According to eq. [5] , the coefficients T are dependent of the 
35 sea bottom geometry and elastic parameter. In the case where 
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these information are not available or are not sufficiently 
accurate, an optimisation similar to the one described in the 
International Patent Application WO96/20417 to determine these 
coefficients numerically- Because the wavefield separation used 
5 here is based on the representation theorem instead of plane 
wave decomposition, the sea floor can be arbitrary irregular in 
the computation of upgoing wavefield described in eq. [4] . 

By applying a scattering series approach, further multiples can 
10 be removed from the upgoing wavefield. 

Though the equation [4] for computing the field {p u , v x u , v y u , v z u } 
removes a significant amount of multiples, including all water 
body reverberations and all receiver peg-leg multiples as 
15 described in FIG. 2A, the process can be extended: Using an 

inverse scattering type operation, as described for example in 
the UK Patent Application GB-A-9710435 , source peg-leg multiples 
(FIG. 2B) can be removed. The first terms of a scattering series 
for the upgoing pressure field are: 

20 

p u (x, CO) = p u (x, CO) + A 1 (CO)?, 17 (x, \, CO) + 

[7] 

A 2 (co)p 2 u (x, \, C0>h . 

where 
25 [8] 

p n u (x s/ £, co) = J dxp u (x, co) J dx'g p (x, co, x')p n _i U (x' , co) . 

Referring now to FIG. 4, the basic steps are described for 
applying the multiple attenuation described above to shot 
30 gathers. Depending on the wavefield from which multiples are to 
be removed, a set of two or more components of the wavefields 
are selected. One of the following combinations must be 
selected: {p, v 2 } , {v x/ v z } , {v y , v z }, {p, v 2 , v x , v y } . Any 
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combinations of the horizontal component can be formed with the 
pressure field or the vertical velocity. The main steps for the 
application are described: 

5 1. Collect mult i -component prestack seismic data. 

2. Select the component to be demultipled and select the other 
components needed for the substraction . For example, to 
demultiple the pressure p, select v z . 

3. Determine the demultiple filter by optimization or by directly 
10 using the equation [5] or both. 

4. Perform the sum in [4] to obtain the demultiple p u . 

An extension of the present invention comprises the step of de- 
composing the measured wavefield into compressional (P-) waves 
15 and shear (S-)wave using, for example, the method described 
above by Amundsen and Re it an. 

For the scope of the invention it is important to note that the 
filter Tof eqs . [4,5] is derived in the x-co domain. Using 
20 standard transformation methods, the filter can be transformed 
into other domains, such as f-k or T-p . 



CLAIMS 



A method of reducing multiples in data recorded during a 
marine seismic survey, comprising the steps of 

(a) providing seismic wavefield data recorded by seismic 
multi -component receivers wherein said data comprise 
measurements related to pressure and velocity components of 
said wavefield; 

(b) choosing a component of the wavefield; 

(c) selecting two or more components of the provided 
wavefield data necessary to attenuate multiples in said 
chosen component; 

(d) filtering one or more of said selected components using 
a spatial filter which combines data components measured at 
different receiver locations and operates indiscriminately 
on P- and S-waves ; and 

(e) summing the selected and filtered selected components 
to generate a reduced multiple variant of said chosen 
component . 

The method of claim 1, wherein the filter includes 
wavefield contributions with zero and non-zero angles of 
incidence at a given receiver location. 

The method of claim 1, further comprising the step of 
decomposing the wavefield into P- and S-waves using a 
separate filter. 



57 .256 



- 15 - 

4. The method of claim 1, wherein the filter is frequency 
dependent . 

5. The method of claim 1, wherein the filter is spatially and 
frequency dependent . 

6. The method of claim 1, wherein the chosen components are 
chosen from a group consisting of pressure p and three 
orthogonal velocity v x , v y , v y . 



10 



7. The method of claim 1, wherein the chosen components are 
chosen from a group consisting of pressure p and three 
orthogonal velocity v x/ v y , v z and the components necessary 
to attenuate multiples in said defined components are 

15 selected according to the following steps: 

if the pressure p is chosen then select the pressure p and 
the vertical velocity component v z to reduce the multiples; 
if the vertical component v z is chosen then select v 2/ the 
pressure p, and the two horizontal velocity components v y 

2 0 and v y ; and 

if one of the two horizontal velocity components v y and v y is 
chosen then select said horizontal component and the 
vertical component v 2 . 

25 8. The method of claim 1, wherein the filter is based on a 

solution of the elastodynamic representation theorem. 

9. The method of claim 1, wherein the filter reconstructs the 
downgoing wavefield from mult i -component measurements of 
30 the total wavefield. 



10. The method of claim 1, further including the step of using 
a correction to adapt for the output of different receiver 
types and sea floor coupling. 

35 



The method of claim 1, wherein the reduced multiple variant 
is further processed using known seismic data processing 
methods . 

The method of claim 1, wherein the wavefield data are 
recorded by receivers located at the sea floor. 

A marine seismic survey, comprising the step of laying out 
a receiver carrier on the sea floor, said receiver carrier 
comprising a plurality of mult i -component receivers able to 
receive at least two independent components of velocity 
and/or pressure measurements; generating acoustic energy to 
travel through earth formations below said receiver 
carrier; using said receivers to measure a mult i -component 
wavefield data of said acoustic energy at receiver 
locations; and the further steps of 

(a) choosing a component of the wavefield; 

(b) selecting two or more components of the wavefield data 
necessary to attenuate multiples in said chosen component; 

(c) filtering one or more of said selected components using 
a spatial filter which combines data components measured at 
different receiver locations and operates indiscriminately 
on P- and S -waves ; and 

(d) summing the selected and filtered selected components 
to generate a reduced multiple variant of said chosen 
component . 
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